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Schrodinger eigenproblems on a discrete interval are further investi- 
gated with special attention given to test cases such as the linear and 
Rosen-Morse potentials. In the former case it is shown that the char- 
acteristic function determining the eigenvalues is a Lommel polyno- 
mial and considerable space is devoted to these objects. For example 
it is shown that the continuum limit of the determinant is obtained 
by a transitional limit of the Lommel polynomials for large order and 
argument. 

Numerical comparisons between discrete approximations and (known) 
continuum values for the ratio of functional determinants with and 
without the potential are made and show good agreement, even for 
small numbers of vertices. 

The zero mode problem is also briefly dealt with. 
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1. Introduction 

In the course of a previous work, [1], investigating the, rather basic, problem 
of the analogue of the Schrodinger equation on a discretised interval (or 'path'), 
Chebyshev polynomials appeared as unperturbed propagators. The use of these 
polynomials in the context of the free equation occurs in earlier works, [2-4]. In 
this communication, I wish to enlarge on the technique employed in [1] of replacing 
the three-term recurrence by a two by two matrix, two-term one, as outlined in the 
classic Atkinson, [5]. There is nothing new in content, but it has a manipulative 
advantage, allowing basic properties and computations to be expressed compactly. 

In addition to developing this formalism (in sections 2 and 5) I present some 
numerical studies of particular potentials comparing the 'exact' continuum results 
with the discrete 'approximations'. I treat the discrete, confined, linear potential 
in some analytical detail and show that the characteristic polynomial is a Lommel 
polynomial. These quantities have not had many physical applications and so i 
present some basic, and not so basic, facts. 

2. The recurrences 

The basic eigenvalue three-term recurrence, (e.g. [6]), 

y(j + 1) + (A - V(j) - 2) y(j) + y(j - 1) = , (1) 

subject to two-point boundary conditions, say Dirichlet, Neumann or Robin, 
y(0) = y(p+l) = 0, D 

y(0) = y(l), y( P ) = y( P +l), N (2) 
Ay{0) = ay{0), Ay{p) = -(3 y{p + 1) , R. 

at the boundary, <9Z, of the discrete interval, X, which is comprised of the (p + 2) 
points j = 0, 1, . . . ,p,p + 1. dX = Up + 1. 

This is a discrete Sturm-Liouville problem. 

V(j) is the potential because (1) can be rewritten as the more familar looking 
Laplacian eigenvalue equation, 2 

[-VA + V(j)]y(j) = Xy(j). (3) 

2 V is the backwards difference operator. There are many discrete approximations to the con- 
tinuum —y"(x) + V(x) y(x) = \ y(x). 



1 



As well as the homogeneous equation (3), the inhomogeneous one 



[-VA + V(j)-\]y(j) = p(j). (4) 

is of interest where p(j) is a source density and A a spectral parameter. 

As intimated in the Introduction, it is helpful to recast these three-term rela- 
tions by two-term matrix ones and I have to use a little space to set this up. 

For example (1) (i.e. (3)), and (4) can be replaced by the matrix recurrences, 

r(j) = C( Xj )T(j-l)^C(j)T(j-l) (5) 

and 

T(j) = C(j)T(j-l) + DR(j-l) (6) 
respectively, with the definitions 

m = ( ?®), ru) = ( , p ®) (7) 



and 

The driving matrix, C, can be split as C(j) = B(j) — \D(j) where, 

/ 1 \ /0 0\ 

B(j) = , D(j) = D = . (9) 

\-l V(j) + 2j \0 l) 

If T x and T 2 are two general solutions of the difference equation, (5), the 
constancy of their Casoratian, 

T 1 (j)JT 2 (j) = det(T 1 ®T 2 ) (10) 

follows immediately, where J is the symplectic metric, 

V-l 

The replacement of a second order difference equation by two, first order ones, 
organised into matrix form, is the simplest example of the more general case where 
an nth order equation is replaced by an n x n matrix equation of the first order. 
Birkhoff , [7] , for example, makes extensive use of this formulation. Milne-Thomson, 
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[8], §5. 3, §12.8, also employs the same technique in the language of continued frac- 
tions. See also [9] and the compact discussion in [10] extended in Agarwal, [11]. 

Taking the discrete interval X as the archetypal domain, the propagation (5) 
can be 'solved' as, 

T(j) = K(j,j')r(j') = K(j,0)T(0), 
where T(0) is the given initial vector and K(j,j') is the forwards propagator, 

k(j,j') = e(j,j') c{j) cu - 1) . ..c(f + 1) , (ii) 

which I concentrate on. 
Equivalent to (5) is 

K(j,j')=5 jjI 1 + C(j)K(j-1,j / ), (12) 

incorporating the initial condition K(j,j) = 1. In particular, for K(j) = K(j, 0) 
one has 

K(j) = 8 j0 l + C(j)K(j - 1) . (13) 

In my previous work, I treated the potential, V(j), as a perturbation. This 
amounts to extracting the factor of 2 — A in C(j) which produces a rearrangement of 
(12) with Chebyshev polynomials as unperturbed propagators. This is formally at- 
tractive, but, if all one is interested in is a numerical answer, then a direct iteration 
of (13) is adequate. Even though it is absolutely equivalent to (1), it is compu- 
tationally advantageous. Furthermore, in order to implement Dirichlet boundary 
conditions, say, it is sufficient to find the bottom right-hand component of K(p + 1), 
i.e. Tr (K(p + 1) D), and the relevant recurrence is 

K D (j) =D5 j0 l + C(j)K d (j - 1) , (14) 

where K D (j) = K(j) D. This is easily programmed. 



3. An explicit example. The linear potential. 

In the continuum, Kirsten and McKane, [12], as a simple illustration of the 
Gel'fand-Yaglom procedure, compute the ratio of Laplacian functional determi- 
nants with and without a linear potential, specifically for the operator, 

L{x) = ~ + V(x), (15) 
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with V(x) = b 3 x. They chose 6 = 1, but I use a general strength. 

As is well known, Airy functions arise for such a potential and Kirsten and 
McKane give, on this basis, 1.085(339648) for the determinant ratio. 

I outline some details later in this section but first I give the discrete version of 
this computation. I proceed to write down the relevant recurrence in scalar form, 
i.e. equation (1) with the potential linear in j, V(j) = Bj, 

y(j + 1) + (A - Bj - 2) y(j) + y(j - 1) = , (16) 

with, for simplicity, Dirichlet boundary conditions at j = and j = p + 1 according 
to (2). 

One might consider this system as describing the motion of a charged particle 
in a uniform electric field, confined to a box, i.e. interval. 

This recurrence is an example of a more general type analysed by Boole, [13], 
using operator methods and, more analytically, by Barnes, [14], but, again I will 
not use their results. Rather, in the particular case of (16), I refer, at the moment, 
to Bleich and Melan, [15], who give a rather detailed treatment of recursions of the 
form, 

y(J + 1) - y(j) + y(j - 1) = , (17) 

i.e. (5) with (8). 

Starting from a solution Ansatz suggested by the iteration of (17), or (5), they 
show that, for a necessarily restricted set of functions, 4>(j), explicit solutions are 
possible. 

In particular a Dirichlet solution for <f>(j) = A + B j is presented in equation 
(59) on p. 149 in [15]. In my notation A = 2 — A and I set their initial point, rj to 0. 
For the moment, I just quote their forwards solution which transcribes (with some 
inessential notational adjustments) to 3 , 

y(j, A) = £ cos ^ ( j - V* ~ 1 ) 3 n 1 [A + B(k/2 + I)} . (18) 

I make further comments on this expression in sections 6 and 7 where I relate 
it to Lommel polynomials. 

To give a flavour of the structure of this solution, I compute y(4, A), which 
corresponds to the terminal value if I choose p to equal 3. 

3 I draw attention to the definition of the product symbol in [15] on p. 141 and footnote on p. 151, 
which differs from the usual one used here. 
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Because of the cosine, k can only be even and therefore only 2 or so that, 



y(4,A) = (2-A + 2S)((2- X + B)(2- X + 3B) -2) 

= -A 3 + A 2 (6 + S x ) - A(10 + S 2 + + 4 + AS X + 2S 2 + S 3 -4B, ^ 

where, for comparison purposes, I have defined v 1 = IB, v 2 = 2B and v 3 = 3B as 
the values of the linear potential Bj at j = 1, 2, 3 and also where 

S 1 ! = + v 2 + v 3 , S 2 = + v x v 3 + v 2 v 3 , S 3 = V X V 2 V 3 . (20) 

I have done this in order to compare with the expression derived in [1] for a 
general potential. The comparison is exact if one takes into account the special 
relation here that VAv 2 = v 1 — 2v 2 + v 3 = for a linear potential. 

The eigenfunction, for a given eigenvalue, say A n , is given by y(j, A n ), (18), 
where n ranges over p values, say 1 < n < p, (setting the terminal point to j = p+1)- 

It is interesting to note that if p is odd, there is always an eigenvalue, halfway 
along, A( p+1 )/ 2 = 2 + B(p + l)/2 for which, as is obvious and can be checked, first 
order perturbation theory is exact. 

Finally, from general theory, the ratio of determinants is given by, 

det(S) y(p+l,0) 



det (B) = rniTZ = »v till (21) 
1 } det(0) y(p+ 1,0)| B=0 ' 

which can be evaluated from (18). 

In order to compare with the continuum linear potential, V = b 3 x, the strength, 
B of the discrete potential has to be rescaled. Over the unit interval, one has 
B = (b/(p + 1)) by dimensions. 

For b = 1, I get a determinant ratio value of 1.08533860, choosing p = 300. 4 
This is not particularly efficient, but comparison with the continuum value, 1.085339648, 
demonstrates the validity of the methods. More interesting, perhaps, is that for just 
one interpolating step, i.e. p = 1, the value is 1.0625, for p = 3 it is 1.07947 while 
for ten, I find 1.08456, a reasonable approximation. 

These results suggests that, if one is prepared to sacrifice a little accuracy, the 
determinant ratio for any potential can be calculated almost by hand. 

Incidentally, for an attractive potential (b < 0), bound states arise and in 
Fig.l I plot the determinant ratio that demonstrates this fact by the zeros. Not so 



4 In 1.4 sees using the wxMaxima CAS on an Athlon ii x 4 machine. 
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interestingly, for a repulsive potential, the ratio shows a rapid, monotonic increase. 
The continuum curve is also shown in the figure. 

Hg.1 . Linear potential determinant 




-7 -6 -5 -4 -3 -2-10 1 

strength b 

For the continuum situation, the Gel'fand-Yaglom function for the potential 
b 3 x is, cf [12], 

y x (x) = ^(Bi(0) M{bx) - Ai(0)Bi(6x)) , (22) 
and yields the ratio of determinants, 

£(Bi(0)Ai(6))-Ai(0)Bi(&)), (23) 

which is also plotted in Fig.l. The figure also shows that the discrete approach is 
good for smallish strengths, b, as expected since it is a perturbation-like method 
yielding a finite polynomial (in b) approximation for the Airy expression, (22). This 
is discussed further in sections 6 and 7. 



4. The Rosen— Morse potential 

Bound states also arise for the Poschl-Teller, hyperbolic potential, 

FW = -^, (24) 
cosh (x) 

confined to a box, as considered by Dunne, [16]. 

From its quantum mechanical origins in molecular physics, the strength, I, is 
usually taken to be integral, but can, in fact, be any real number, in which case it 
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might be preferable to use the variable 1 = 1 + 1/2. I will still employ / but refer to 
(24) as the Rosen-Morse potential (e.g. Kleinert, [17], Fliigge, [18]). 

Two particular solutions of the equation, L(x)ip(x) = 0, 5 are the Legendre 
functions, P ; (tanhx) and <5 ; (tanhx), constituting a fundamental set. I choose the 
interval —1/2 < x < 1/2 as domain so that I can give the Gel'fand-Yaglom function 
as written down by Dunne, 

ifj GY ( x ) =P l (- tanh 1 /2) Q t (tanh x) - P l (tanh 1 /2) Q t (- tanh x) 

so that the ratio of determinants, with and without V, is 



det c (I) = P^-t) Q t (t) - P t (t) Q t (-t) (25) 

where t = tanh 1/2 0.46211715726. Properties of the Legendre functions check 
that det c (l) is symmetrical in I, as must be, and one need compute it only for 
positive Z, i.e. I > —1/2. 

In terms of functions of the first kind, I find, 

te c {l) = -^(P l {-t)*-P l (t)*) (26) 

which is easily calculated using Gauss' hypergeometric form, say, 

P l (t) = 2 F 1 (-l,l + l;l;(l-t)/2). 

Near integral /, the numerics become uncertain. But for / actually integral, the 
ratio of determinants, (25), can be expressed purely in terms of the polynomials, 

Pi, 



det c (I) = (-l) l P l - 2 ^ P m-l P n -m) , ^ Z , 

^ m=l ' 



(27) 



which is a polynomial in tanh 1/2 and provides a check of the numbers. As an 
exercise, one could use these values as the basis of an interpolation. 

Turning now to the discrete calculation, the formula I use is that developed in 
[1]. For p = 3, an interval with three interior vertices and step size h = 1/4, the 
determinant ratio is, 

1 



det D (l) = l + -(SS, + 2S 2 + S 3 + v 2 ) , (28) 



5 For simplicity, I set the mass equal to zero. 
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where the S i are given by (20) with Va = V(j). (Refer to (3)). This is to be 
compared, numerically, with (26). 

On this side of the calculation, I use the interval to 1 so that the discrete 
potential corresponding to (24) is given by 

= h*l(l + l) 
cosh 2 (/ij-l/2) 

and the evaluation of (28) is relatively elementary. The results are presented graph- 
ically in fig. 2, where I have also displayed the 5 vertex values and the continuum, 
'exact' curve, based on (26). 

The figure shows again that the discrete approach is good for smallish strengths, 
I. It yields a finite polynomial (in I, even in 7) approximation for the Legendre 
expression, (25). The exact curve oscillates for ever and, while a finite polynomial 
cannot do so, it follows fairly well, for a decent range of Note also that, for a 
repulsive potential, the ratio is bigger than one. 

Fig. 2. Rosen-Morse determinant 



: 1 : 

continuum 




2 4 6 8 10 12 
strength I 

The zeros of the curves occur when bound states appear as the well deepens. 

5. Analytical extensions 

So far it has been assumed that label j takes just integral values. However in 
classical finite difference theory, e.g. Boole, [13], this is only a particular case of a 
more general situation where one deals with functions y(x) of a variable, x (real or 
complex) satisfying a difference equation, say, 

y(x + l)-<f>(x)y(x) + y(x-l) = 0, (29) 
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as a relevant, typical example. In accordance with my previous formalism, I write 
this, with, therefore, some repetition from section 2, as, 



T(x) = C(x)T(x-l) 



(30) 



with 




(31) 



and 




(32) 



Assuming some starting point, the general solution of (30) is 



T(x) = C(x)C(x - 1) . . .C{x + l)T(x ) 
= K(x,x )T(x ) 



(33) 



where T(x ) can be considered as an arbitrary constant vector, say O. 

For a particular x , the domain of x is determined, but one can consider x to 
change, and with it, the domain of x, (see e.g. Boole, [13], p. 102, Levy and Lessman, 
[19], pp. 89+) (always, however, with x — x Q integral). This can be accommodated 
by making the constant, O, a function, 



where the w i are arbitrary unit periodics, Bohmer, [20], [15], §26, [19], [8], §11.1, 
(e.g. w^Xq) = w^x)) so that the general solution, (33), is, 




T(x) 



C(x)C(x - 1) . . . C(x + 1) fl(x ) 

K(x,x )tt(x ) 

zu a (x)T a (x) + w p {x)Tp{x) , 



(34) 



where, 




K(x, x )ct , T p(x) = K(x,x )f3 , 



(35) 



with the standard basic vectors, 




(36) 



For example, Dirichlet initial conditions require w a (x ) = 0. 
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The solution (34) is just the general solution to the second order difference 
equation (29), a standard construct, e.g. Levy and Baggott, [21], Birkhoff, [7]. 

The basis, T Q and is a canonical one {e.g. Lakshmikantham and Trigiante, 
[10], Lemma 2.1.1). 

T p coincides with Atkinson's standard solution, y n (X), [5], equns. (4.1.4), 
(4.1.5) and T Q with the standard solution, z n (\), equns. (4.2.6), (4.2.8). 

As a non-trivial example, the solution (18) to (17), corresponds to the Dirichlet 
solution Yp{j) because y(0,X) = and y(l,A) = 1. In the following section, I 
enlarge on this solution. 

The constancy of the Casoratian, (10), now becomes the statement that it is a 
unit periodic, 

T 1 (x)JT 2 (x) = w(x) , 

where T l and T 2 are two general solutions of the difference equation, (30) each of 
the form (34), say. 

6. Bessel functions 

I return to a consideration of the confined linear potential of section 3, which 
forms an interesting, if somewhat specific, case but with connections to some special 
functions. 

The analysis of the second order equation by Barnes and by Bleich and Melan, 
proceeds from the classic form of the recurrence, (29), so I present their results 
accordingly, despite my preference for the matrix form. 

The equation in question is (cf (16)), 6 

y(x + 1) - (Bx + A) y(x) + y(x - 1) =0, (37) 

which, treated using the textbook Laplace method, 7 yields two independent solu- 
tions as contour integrals coinciding, to a simple factor, with Bessel functions, 

yW(x) = (-1)*+*/* J_ (x+A/B) (2/5) = J. {x+A/B) ( - 2/B) 
y^(x) = J {x+A/B) (2/B) 1 



This is recognised as a recurrence relation for Bessel functions but it is not necessary to do so 
at this time. 

A simple discussion is given in the elegant little book by Miller, [9] . 
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and the general solution is 



y(x) = w x {x) y^\x) + w 2 (x) y^(x) . (39) 

If £ = x + A/B is an integer then these two solutions are equivalent and, in 
the usual fashion, a second solution is obtained by a limiting process, £ — > n, on a 
Bessel function of the second kind, due to Weber and Schlafli, 

Y^z) = cot frJ^z) - — L- J_^z) , (40) 
s s sin^7r s 

which satisfies the difference equation (29) as it is of the form (39). 

As an application, return to the conventional case when x = j e Z with the 
previously used discrete interval, X, as domain and attempt to recover earlier results, 
such as (18). Thus, apply D conditions at j = and j = p + 1. This produces the 
non-local Casorati determinant condition, 



yW(0) y( 2 )(0) 

y (D(p+l) ^(p + l) 



, (41) 



which is just the vanishing of the Dirichlet solution, 

at the terminal point, x = j = p + 1. 

Inserting the explicit functions, (38) gives the quantity, 8 

W(u,p) = J_J p+1+u -{-iy +1 J_ {p+1+v) J„ = 0, P GZ, (42) 

where u = A/B and the arguments of the Bessel functions are all 2/B = z . I 
have dropped inessential, common overall constants and leave the dependence of 
W, which is hereby defined, on z implied, until later. 

The characteristic equation, (41), in A (remember, A = 2 — A) provides the 
eigenvalues through its roots, and is easily computable since, in this case, two 
independent solutions are known. 

The aim, then, is to show that the Casoratian, (41), is proportional to a finite 
polynomial in v of order p, cf (18). This can be done via the series form of the 
Bessel functions in (42). Of course, it is to be anticipated that the Casoratian 



The relation is y(p + 1,A) oc W{v,p). 
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vanishes when v is an integer for then the two Bessel functions are proportional. 
To allow for this, it is convenient to consider the normalised Casoratian, 9 

— ._W(u,p) 
W(u,p) = 



W(u t 0) ' 



which removes the relevant factor. 
A standard formula is, 10 



W(is, 0) = J_ U J 1+U + J-i-ty J v = sin un . (43) 

7T 

I give the calculation and refer to comments in sections 7 and 8. The compo- 
sition relation is used to give the series, 

^ oo^ (_ 1 )n g -(2n+p+l)( p + n + 2 ) 



n=0 



(-l^B-^-P-^in-p), 

J -p-l-v — /^l 



— ' n\T(n + 1 + u)T(n — p — v) 



n=0 



In the combination, (42), most of the terms cancel, leaving the finite sum, 

W<y p) = y( P + 1 - 2n ) ( - 1)ng " (P+1 " 2n) (44) 



n=l 



Part of the summation range gives zero because (^ p+ ]_~ 2n j = , if n < (p+ l)/2. 

I give the details that turn (44) into (18), up to a factor. The Gamma functions 
combine as, 

1 (u + p-n + 2)...(u + n-l) . 

T(-u -n+l)T(p-n + v + 2) tt 

(A + B(p - n + 2)) . . . (A + B(n - 1)) 

= kB*"-p-* SmU7V 

(45) 

and the normalised Casoratian is, 
W(A, = £ ( P + 1 _~ 2n )(-l) n+1 {A + B(p-n + 2))...(A + B(n - 1)) , 

n=l ^ ' 

(46) 



9 An alternative procedure is given in the next section. 
10 p = means there ar no (internal) vertices. 
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The number of brackets being (p — 2s). 

I have changed arguments to allow for the z dependence. Thus, W(A, B,p) = 
W(v,p)(z). 

The change of summation variable to, 

k = 2(p + 1 -n) , 

then yields agreement with (18) (if j is chosen to be the end boundary point, p + 1) 
since, 

/p+l-2ra\ = +1 _ n f n-1 \ 
V -n J 1 ' \2n-p-2)' 

Bleich and Melan's solution, (18) has thus been obtained by a more direct and 
particular route. 

Finally the trivial transformation of summation variable s = p + 1 — n converts 
(46) into, the simpler-looking form, 

<(p+l)/2 / _ x 

W{A,B,p) = {-iy V S )(-l) s (A + B(p-s))...(A + B(s + l)), 

s=0 V s / 

(47) 

where the number of brackets is equal to (p — 2s). 
The case p = 1, 

Jv+2( Z ) J-v( Z ) ~ J-(v+2)( Z ) Jv( Z ) = " \~ ) 

7T Z 

is given in Gray and Matthews, [22], p. 241, Ex.1. Other examples are, 11 
T (A T M4- T (A T (A ^ ~ ^ + ^ + 2)) smun 

Jv+±( z ) J-A z ) - J-(v+4)( z ) JA Z ) = " — z 4 • 

(48) 

One could regard the solution (18) as an alternative source of these results 
which actually are most easily obtained by direct iteration of (37) to give T 2 (j) of 
(35) as described in the next section. 

I now place these evaluations in a more general and historic setting. 



11 For odd p, the vanishing of the overall factor (y + (p + l)/2)) corresponds to a cancellation 
between the two terms on the left-hand side which does not require using J n = J_ n . 
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7. Lommel polynomials 



There is no obligation to identify solutions of the recurrence equation, (37), 
specifically with Bessel functions but making the identification saves a lot of effort 
due to an extensive literature. For example, the algebra leading to (44), and the 
further extraction of the factor smvn, is well known. What I have termed the 
normalised Casoratian is a Lommel polynomial (Lommel, [23], Nielsen [24], Graf 
and Gubler, [25], Watson, [26]) as I now show. 

Lommel defined the polynomials, R, from the iteration of the particular Bessel 
recurrence, 

VW = 7^)- f mW ( 49 ) 

specifically for F^(z) = J^(z) (not the only choice). 

Iteration yields the formal solution in terms of two initial values as the reduction 
formula, 12 

i.e. 

F^z) = R^-^(z) F c _ n (z) - Rt-w-^z) F^^iz) , (50) 

where, to confirm the connection with my previous notation in (37), z = 2/B 
£ = x + A/B and y(x,A,B) = F x+A / B (2/B). I have included the parameters A 
and B in the y(x). Either £ or x can be used as the current, discrete coordinate. 
Writing (50) in terms of the initial point £ , i.e. eliminating n = £ — £ — 1, 

F € (z) = fl«o,€-€o-i( z ) F€o+1 ( z )_ jR €o+i,€-€o-2 (z)F€o ( z ). (51) 

Comparison with (35) shows that i? ? o. ?-?o- 1 ( /2 ) corresponds to the /3, or 
Dirichlet, solution, and —Rt° +1 't~t°~ 2 (z) to the a one, as functions of £. Checking 
this requires use of Graf's relation, [25], 

p^-n-l ■p£v—n,n—l 

which implies R v ~ x = and R v ~ 2 = -1 e.g. [25], p.103. 

Lommel deduced the actual form of the polynomials, R^' n (z), by an inductive 
procedure repeated by Graf and Gubler, [25] vol.1 chap.l §2. A simpler method, 

12 I use Nielsen's notation, [24], which is the same as that of Graf and Gubler, [25], who write £ . 
It differs from Lommel's and Watson's only in the ordering of the indices and a shift of unity in 
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e.g. Nielsen [24], is that give in section 5, since easy algebra leads to Lommel's 
equation, 



J-v(z) J„ +P+ i(z) + ("1) P ■/„(*) J-u-p-M = ~ 2 -^P R"' p (z) , (52) 



and the left-hand side can be evaluated from the known series form of J v {z). The 
details in section 5 are, more or less, just a rewriting of Nielsen's algebra and I also 
give the standard result from this old literature, [23], [24], 

identical to (47) . 

A comparison of (52) with (42) gives the connection, 

W(u, p) = - 2SH17n/ R"*(z) , z = 2/5, 
nz 

and the normalised Casoratian is then identical to the Lommel polynomial, 

W(v,p) = R u ' p (z) , 

for non-integral v. This extends to integral v by continuity since R n,p is uniquely 
defined (cf Watson, [26] §9.61). 

The classic result, (43), is a consequence of the general (52) for p = together 
with the explicit form of the R u,p , although it can be established independently. 

In order to regain the notation of (37), v must be set equal to A/B, and z to 
2/B and (53) yields the most explicit formula (47), as in section 5. (See below.) 

Actually, the passage to integral v is best approached using the second kind 
functions, Y v . Simple algebra using (40) shows that (52) is replaced by the more 
elegant, 



— R v *{z) = Y v {z) J v+p+1 {z) - J u (z) Y u+p+1 (z) = 

77 Z 



Y v Y u+p+ i 



(54) 

also due to Lommel, [24]. This will play a role later in the continuum limit. 
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8. Lommel polynomials as orthogonal polynomials 



First a few words regarding variables are relevant here. The variables that 
provide the standard form of the Lommel polynomials, i.e. equation (53), are z 
and v related to the 'physical' variables A = 2 — A and B where A is the eigenvalue 
and B the strength of the (linear) potential, by v = A/B and z = 2/B, as stated 
earlier, so that 

2-A z , n ^ 
"= — =2 (2 - A) 
In terms of A and B separately, (47) provides the expression for Lommel poly- 
nomials. It shows that they are polynomials in both A and B or extracting factors 
of A or of B, polynomials in B j A or A/B = v. Here I am interested in polynomials 
in A = 2 — A, or equivalently for fixed B, in v. 

At this point, it is convenient to list a few Lommel polynomials, 13 



R(v,l) = (u+l) 



2^ 2 



2\ 3 , J2 



(55) 



R(u,2) = (u + 2)(u+l)^- 

R(u, 3) = (u + 3)(u + 2){y + 1) (JJ - 2{v + 2) ^ 

i2(i/,4) = (i/ + 4)(i/ + 3)(i/ + 2)(i/ + l)0^ -3(^ + 3)^ + 2)0^ +1 
together with the associated recursion (due to Lommel, [23]) [24], p. 54, [25] p. 102, 

R(v,p + 2)= 2 fr +I/ ) R(y,p+1)-R(u,p). (56) 

z 

By classic, general theory, the y(j,X) of (18), or equivalently the W of (46) 
(with A = 2 — A and p + 1 reset to j), i.e. the Lommel polynomials, R, form a finite 
set of explicit orthogonal polynomials (in A, or v) concentrated on the eigenvalues, 
X n or u n , e.g. Atkinson, [5]. 

For p vertices, orthogonality reads 

v v 
£ R(v, R{v\ j~l) = l + Yl j' 1 ) ^ 3-l) = N(*) S »w ( 57 ) 

J=l 3=2 



13 There is a misprint on p. 102 in [25]. I have temporarily changed notations so that R(v,p) = 
R u >p(z). 
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where v and v' are eigenvalues determined by the vanishing of R(u,p). 

It is interesting to confirm this relation. As a numerical example, take p = 2. 
The eigenvalue equation is then 

(z/ + l)(z/ + 2)-^ 2 = (58) 

and (57) reads 

1+ + + = 0, v±v' (59) 

Z 

and 2 

l + fL(i/ + l) 2 = tf(i/). (60) 

Using [y + 1) as the variable, the product of eigenvalues is, from (58), —z 2 /4 and 
(59) is readily verified. (60) fixes the normalisation for orthonormality. 
More complicated is the case of p = 3. Orthogonality is now 

1 + %{y + + 1) + [^{y + 2){y + 1) - 1^1 f-U' + 2)(^ + 1) - = (61) 



From (55) the eigenvalue equation R(u, 3) = has a solution v = —2, as 
mentioned before, and, for simplicity, I set one eigenvalue, v' = —2 so that (61) 
reads, 

1 - \{y + 1) + 1 - \ {v + 2){y + 1) = 

which reduces to 

(^ + 3)^ + 1)^-1 = 0, (62) 

z 

i.e. the equation determining the other two roots of R(u, 3) = and the confirma- 
tion is complete. A similar mechanism works for any odd p. 

Showing that (61) holds for solutions of (62) is more complicated. 



9. The continuum limit 



In this section I show how to derive the continuum determinant ratio, (23), 
from the above discrete expressions. 

I denote by h the path step, h = l/(p + 1), and seek to let h — > 0, particularly 
in (54). As noted earlier, the strength, B, of the discrete linear potential is related 
to 6 3 , the strength of the continuum linear potential by B = (b/(p + l)) 3 = (bh) 3 . 
Then the argument, z = 2/(bh) 3 — > oo. 
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Since I am interested only in the determinant, I can immediately set A = 
so that v = z. Therefore as h — > the leading divergence of the the orders of the 
Bessel functions in (54) is just that of their arguments and I can employ Nicholson's 
asymptotic relations. Starting from 

2 R**{z) = Y z {z) J z+p+1 (z) - J z {z) Y z+p+1 (z) (63) 



7T Z 



I need 



n+„w~3- 1/6 (f) iMO+j-vato] 



and 



= © 1/2 3-1/6 (!) 1/3 AiW 



where Ai and Bi are Airy functions and 

H© 1/2Q3/2 = r 3/2 - 

The parameter a is either or (p + 1) = 1/h and 2/z = (hb) 3 from which it is 
seen that the arguments, x, of the Airy functions are either or b. 
The multiplying factor reads 

ox 1/2 /A l/3 /l\ 1 /3/ 9 xl/3 / 9 \l/3 

<M 3-i/e/ £ A =3 i/2 3 -i/6(M /2 



and there are two of these. 

Hence we have the transitional limiting behaviour of the Lommel polynomials 

R z ' p (z) ~ 7T ( - ) (Bi (0) Ai (b) - Ai (0) Bi (&)) 

V 2 / (64) 

~ \ t (Bi (0) Ai (6) - Ai (0) Bi (&)) . 
ft o 

The form (64) contains within it Cauchy's limit, 

r(i/3) 



22/331/6^1/3^1/3 
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and also 




3 1 / 3 r(l/3) 
22/370,1/3 ' 



Numerically, choosing p = 9, a not too big a number, i.e. h = 1/10, the exact 
evaluation of #(2000, 9, 2000) from the series is 10.84393086, for b = 1, whereas the 
limiting approximation, (64), gives 10.85339648. (These are just the same as the 
numbers in section (3).) 

As I have not been too careful with normalisation, this quantity, (64), will be 
proportional to the functional determinant. This does not matter as I only want 
the ratio of the determinants, with and without b. As b — > 0, it is trivial to show 
that the right-hand side of (64) is just l/(hb) and the desired ratio agrees with 
that, (23), obtained in the continuum directly by Kirsten and McKane, [12]. 

10. The discrete zero mode problem 

If one of the eigenvalues, say A , vanishes, then the ordinary determinant is 
zero. It is then customary to calculate the product of the remaining non-zero 
eigenvalues and refer to this as the determinant. Kirsten and McKane, [12], give a 
modification of the GePfand Yaglom procedure that covers this case. Their formula 
just involves the zero mode eigenfunction. In this section I give the discrete version. 

The result drops out of the discrete analogue of Green's theorem, or Lagrange's 
identity, 14 , or the Christoffel-Darboux identity, e.g. Atkinson, [5] §4.2. For speed, 
I copy the theorem from [5], for the recursion (1) 



where A and \i are any two numbers. Now choose \i to be an eigenvalue, A i? and set 
the upper limit k = p, the number of vertices. Then y(p + 1, A) is the characteristic 
polynomial and so y(p + 1, y) = y(p + 1, AJ = 0. Equation (65) is rearranged into 




y(k + l,fj) 
y(k,fj) 





(66) 



Setting A to Aj gives 



1 



v 



MAJ I y{\)) 
Ay(p, \) 



ri( A - A n) 



yip, K) 



(67) 



This is the method employed in [12]. 
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where I have again used the Dirichlet condition y(p + 1,AJ = 0. Choosing A^ = 
A = 0, the left-hand side is, up to a sign, the determinant, omitting the zero mode. 
This result should be compared with the relevant part of equation (14) in [12]. 

As a non-trivial example I again turn to the linear potential case when the 
characteristic polynomial is a Lommel polynomial. Specifically I choose 3 vertices 
(p = 3) and select the v = — 2 root to give the zero mode. This means that the 
strength, B, has to equal —1 since v = (2 — X)/B. The remaining two roots are then 
easily found as ±\^3 so that their determinant is —3. Computation of the relevant 
Lommel polynomial, R(— 2, j — 1)(— 2), gives the zero mode as (1, 1, —1), which is 
consistent with (67). 

Incidentally, by dividing (66) and (67) one obtains Legendre's function and 
Legendre's polynomial interpolation of / can therefore be written as, 



in terms of a set of p orthogonal polynomials. (To repeat, the eigenvalues, A n , are 
the zeros of the polynomial, y(p+ 1, A).) An historic and practically important case 
corresponds to the free equation, V = 0, which yields Chebyshev polynomials. See 
[27] chap.6. 

Equation (68) is related to a discrete Kramers sampling theorem, cf Annaby, 
[28]. If /(A) is a linear combination of the y(j,X), (1 < j < p), then F(X) = /(A) 
and /(A) can be reconstructed from its samples, which is, perhaps, not surprising. 

This can be more abstractly expressed as the vector space expansions, 



(68) being regained on multiplying by (A |, the function label, y, being suppressed. 
The condition on / corresponds to the bandwidth limitation of the original Whittaker- 
Shannon sampling theorem. 

11. Comments 

It might be said that the Bessel functions are only intermediaries in the cal- 
culation of the Lommel polynomials, which could be taken as the primary objects. 
Indeed, Bessel functions can be obtained as limits of Lommel polynomials, as shown 




(68) 
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by Hurwitz, [29]. (Incidentally, precisely this limit was discussed by Bleich and 
Melan, [15], who make no reference to any earlier work.) 

If this attitude is adopted, the form of the Lommel polynomials, e.g. (53), must 
be found without explicit use of Bessel functions, but using only the single recursion 
(37). Such is Lommel's original computation, [23], but he assumes the result and 
justifies it by induction which is somewhat synthetic. Bleich and Melan's method, 
[15], is interesting, but a little complicated, and again involves an assumption re- 
garding the solution's general form. Graf and Gubler, [25] vol.2, however, use the 
recursion, expressed as a continued fraction, and obtain the Lommel polynomials, 
which they refer to as Schlafli functions, directly. They list a few examples and also 
give the factored form (47). 15 

It is possible to derive any particular Lommel polynomial by direct iteration 
of the recursion (49), or (37). The reduction formula, [24], §26, 

R">P(z) = cos V — + - ^(u + s + 1) sin ^^tt R»> s 
2 z s=o 2 

can also be iterated. 
11. Conclusion 

It has been shown that, apart from providing quite accurate approximations for 
smallish numbers of lattice points to a continuum situation, discrete techniques have 
an intrinsic interest and form a useful intro to the theory of orthogonal polynomials. 
In particular, I have highlighted the case of Lommel polynomials which seem to 
have been ignored. Even their general orthogonality properties do not seem to be 
available explicitly. 
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